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The differential Adhesion Hypothesis (DAH) is a theory of the organization of 
cells within a tissue. In this study we introduce a stochastic model supporting 
the DAH, that can be seen as a continuous version of a discrete model of 
Graner and Glazier. Our approach is based on the mathematical framework 
[ of Gibbsian marked point processes. We provide a Markov chain Monte Carlo 

I algorithm that can reproduce classical biological patterns, and we propose an 
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estimation procedure for a parameter that quantifies the strength of adhesion 



O ' between cells. This procedure is tested through simulations. 
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1 Introduction 



The development and the maintenance of multi-cellular organisms are driven by 
permanent rearrangements of cell shapes and positions. Such rearrangements are a 
key step particularly towards reconstruction of functional organs [1^. In vitro ex- 
periments such as Holtfreter's experiments for the pronephros and the famous 
example of an adult living organism Hydra are representative example of spec- 
tacular spontaneous sorting. 

Among the pioneering works, Steinberg used the ability of cells to self-organize 
in coherent structures to conduct a series of experimental studies that characterized 
cell adhesion as a major actor of cell sorting mEHElIZI- Following his experiments, 
Steinberg suggested that the interaction between two cells involve an adhesion sur- 
face energy which varies according to the cell type. He interpreted cell sorting via the 
Differential Adhesion Hypothesis (DAH), which states that cells can explore various 
configurations and finally arrive to the lowest-energy configuration. During the past 
decades, the DAH has been experimentally tested in various situations such as gas- 
trulation [Sj, cell shaping jU], control of pattern formation |10j and the engulfment 
of a tissue by another one. Recently, some of these experiments have been improved 
to evaluate the DAH rigorously |TT| . 

Many mathematical models have been previously developed for the DAH. Most 
of these models are dedicated to the simulation of physical processes in agreement 
with the DAH: model simulations tend to minimize an Energy functional, called 
the Hamiltonian, supporting the DAH. Four groups of models can be distinguished 
according to the geometry of the tissues. First, cell-lattice models own the particular- 
ity that each cell is geometrically described by a common shape, generally a regular 
polygon (square, hexagon, etc ...) (see [12] for example). The second class of models 
has been called centric models. This class provides more realistic cell geometries by 
using tesselations (Dirichlet tesselation for example) to define cell boundaries from 
a point pattern where points characterize cell centers J^l- The third class of models 
is called the vertex models. These models are dual to the centric models jEl 1^ • 
Thr fourth class of models, called sub-cellular lattice models, has been developed as 
a trade-off between the simulation speed of cell-lattice models and the geometrical 
flexibility of the centric models. These models have been introduced by Graner and 
Glazier JH], and they are also referred to as the discrete GG model afterwards. 

In this study, a continuous model for cell sorting derived from the discrete GG 
model is presented. In this context, the geometry of cells is actually similar to the 
centric models: assuming that cell centers are known, the cells are approximated 
by Dirichlet cells. In agreement with the DAH, the new model is also based on an 
Hamiltonian inspired from the GG model Hamiltonian (see Section 12)). Following 



this approach, the model falls into a well-defined mathematical framework: Gibbsian 
marked point processes ^Zj. This mathematical background will allow us to better 
control the simulation procedure for generating cell sorting patterns. 

Furthermore, recent developments in molecular biology emphasize the fact that 
cell-cell interactions play a major role in tumorigenesis ^Hl- The nature of the 
interactions may actually reflect the initiation of a cancer. In addition, the invasive 
nature of a tumor is directly linked to the modification of the strength of cell-cell 
interactions In this perspective, an important challenge is to quantify the 

strength of the adhesion between cells. Another goal of this study is therefore to 
provide an inference procedure for the parameter that governs the strength of cell-cell 
adhesion. 

The article is structured as follows. In Section |21 the continuous model is in- 
troduced. Section El presents an estimator of the adhesion strength parameter, and 
gives some mathematical properties of the model. In Section HI results concerning 
simulations of cell sorting patterns and the performance of the statistical estimator 
are presented. 

2 A continuous cell model for DAH 

This section introduces a new continuous model for differential adhesion. As in 
the previous approaches, the model is based on an Energy functional that describes 
cell-cell interactions. In the continuous model, the tissue is still described by a cen- 
tric model where the points correspond to the locations of the cell centroids, and 
the marks correspond to the cell types. Honda's studies showed that the geome- 
try of Dirichlet cells was in agreement with biological tissues for which the spatial 
coordinates of the nuclei were extracted thanks to molecular markers |2(H l2T] . 

The GG model. Before describing our model, we start by giving a brief account 
on the GG model ^E]- In the GG model, a cell was not defined as a simple unit, 
but instead consisted of several pixels. The pixels could belong to three types: high 
surface energy cells, low surface energy cells or medium cells, which were coded as 
1,2 and -1 respectively. According to the DAH, the Hamiltonian Hqq was defined 
as an extension of the Potts model as follows 

Hgg= ^ i^^^d)^ ^i'^i'd')) (l - + ^ - ^rH)'r(yl,(,)), 

(1) 

where were the pixel spatial coordinates, aij represented the cell to which the 
pixel belonged, r^dij) denoted the type of the cell aij, and the function J char- 



acterized the interaction intensity between two cell types {6 denoted the Kronecker 
symbol). In particular, the term ^1 — ^o-i j,o--, ^, j indicated that the interaction be- 
tween two pixels within the same cell was zero. Shape constraints were modeled 
by the second term where A corresponded to an elasticity coefficient, a{a) was the 
cell area and v4^(o-) was a target area that depended on the cell type. The function 
r denoted the Heavyside function. It was included in the formula so that medium 
cells (coding -1) were not subject to the shape constraint. 

The continuous model. Here we denote by (xi) {i = 1, . . . ,n) the n cell centers. 
The Dirichlet cell of Xi is denoted by Dir(xi), and is defined as the set of points 
which are closer to Xi than to any other cell nucleus. In addition we write |Dir(xi)| 
to denote the area of the cell Dir(xj). A continuous version of Equation^ can be 
constructed as follows. In the GG model, a cell a is in neighborhood of a cell a' as 
soon as a single pixel of a is adjacent to a pixel from a'. With this in mind, the GG 
model's Hamiltonian can be rewritten as 

Hgg = 5^ k n a'l J (r(a), r(a')) + A 5^(a(fT) - A,^,^fT{A,^^)) 

where |(T fl cr'| is the number of connected pixels between a and cr', which can be 
identified as the Euclidian length of the interaction surface between the two cells a 
and a' . Identifying cells to their centers, |crncr'| can be approximated as |Dir(a;jna;j)|, 
where |Dir(xinxj) | is the length of the Dirichlet edge shared by Xi and Xj. In addition 
a cell area in our model matches with the area of a Dirichlet cell, which means that 
a{a) corresponds to |Dir(xj)|. Defining a (marked) cell configuration as 

= {(xi,ri),...,(x„,r„)}, (2) 

where the ( the cell centers, and the (r,) are the corresponding cell types, we 

define here an equivalent Hamiltonian function 

H{^) = |Dir(a:. n x,)| J(r„ r,) + A 5^(|Dir(x,)| - A,^fT{A^^). (3) 

The symbol i ^ j, means that the cells Xi and Xj share a common edge in the 
Dirichlet tiling. As in the GG model, this Hamiltonian can be viewed as the sum of 
two terms. The ffist term corresponds to pair potentials and controls the adhesion 
forces between contiguous cells. The second term corresponds to singleton potentials, 
and is analogous to the one introduced in the GG model. It controls the shape of 
the cells as well as their density. 

An additional parameter, denoted 6, is also included in the model. It is called 
the adhesion strength, and allows us to define the Energy functional of our model as 



6H[(p). Because A is a free parameter, we see that 6 actually controls the relative 
intensity of pair interactions (adhesion between cells). This parameter is of partic- 
ular interest because an important benefit of the continuous approach is to allow 
consistent statistical estimation procedures for this parameter. 

Surface tensions. Biological tissue configurations are often interpreted in terms 
of surface tension parameters. For instance, checkerboard patterns are usually as- 
sociated with negative surface tensions, whereas cell sorting patterns are associated 
with positive surface tensions. When two distinct cell types are considered, the 
surface tension between cells with the distinct types can be defined as 



where ri, T2 are the cell types. In [22], surface tensions were denoted as 7/^. In the 
next section, the interaction parameter will not correspond to those used in the GG 
model exactly. However, they will be fixed so that the surface tensions are of the 
same order. 

3 Model simulation and parameter estimation 

There are two important benefits of assuming a continuous model for cell sorting. 
Following this approach we will be able to 1) provide mathematical conditions for 
warranting the convergence of the model simulation algorithm (such controls were 
missing from the original discrete approach), 2) propose statistical procedures for the 
inference of the interaction strength. To reach these objectives, we shall integrate 
our model in the theory of Gibbsian marked point processes which provides a general 
framework for simulation and parameter estimation (see PT| 123 ) ■ 

Simulation. Simulation from the continuous model can be performed according 
to the Metropolis-Hastings algorithm. At each iteration, the algorithm randomly 
chooses between two operations: either inserting or deleting a cell within a well- 
delineated region. Insertion and deletion of a cell in the configuration has been 
implemented using local changes as described in |21j and [25.. In both cases, the 
algorithm computes the Energy variation after one operation is performed. Then 
it may accept the operation with probability p = min(l, exp(— ^^A//)) otherwise it is 
rejected. Convergence results for this continuous state space Markov chain will be 
established afterwards. 




Pseudo-likelihood inference. In this section we propose an inference procedure 
for estimating the parameter 6. Estimating 9 actually provides information about 
the strength of cell-cell interactions and adhesion. Here we resort to a classical 
approximation in statistics: the pseudo-likelihood method, first introduced by Besag 
in the context of the analysis of dirty pictures [211 • For any configuration ip, the 
pseudo-likelihood is defined as the product over all elements of ip of the following 
conditional probabilities 

PL(^^)= J] Prob((x„r,)|¥P,^) 

In this formula, the conditional probability of observing (xj, Xj) at Xi given the con- 
figuration outside be described as 

where M corresponds to the set of the possible cell types (or marks), and where 
H^{xi,Ti) represents the contribution of the cell Xi in the expression of the Hamil- 
tonian H{{p), i.e.. 

We also consider the logarithm of the pseudo-likelihood given by 

L(^) = - ( ^H^{xi, Ti) + log exp{-eH^{xi, m)) J . 

Maximizing Equation ^ provides an estimator of 6, namely 

9{ip) = a.TgmaxgL{ip,9), 
which can be achieved using standard techniques. 

Gibbsian marked point processes. We now recall some basic results about 
Gibbsian marked point processes. Gibbsian models, according to the statistical 
physics terminology, have been introduced and largely studied in j22j or [2HI. The 
Gibbsian category is analogous to the Markov point processes introduced in spatial 
statistics by Ripley and Kelly [SH] and reviewed in great details by Van Lieshout 

Here, we restrict ourselves to marked point processes that have a density / with 
respect to the standard Poisson process. A realization of such a process is called a 
configuration and is denoted as (p. When cp has exactly n points, we can write 



(4) 



ip = {{Xi,Ti), {Xn,Tn)}, 



as in Equation [21 A cell-mark couple (xj,rj) is called a point. According to the 
Hammersley- Cliff ord-Ripley-Kelly Theorem (see jSHj), the density of the process 
conditional to the number of points is of the following form 



where H is the Hamiltonian of the system, Z is the partition function. Using this 
formalism, the convergence of the Markov chain generated by the above algorithm 
can be studied. We focus on the Harris-Recurrence property, which means that the 
convergence of the chain is ensured for any initial configuration with a non-zero 
probability. The following result can be stated. 

Proposition 1 Let us consider a Gibbsian marked point process as defined in Equa- 
tion\^ and 



where J charaterizes the interaction intensity and A is an elasticity parameter. 

Assume that J is bounded (\J\ < Kj) and that two neighboring points are at a 
distance between d and D (0 < d < D < +oo^. Then the Markov Chain generated 
by the simulation algorithm is Harris-Recurrent. 

The proof of proposition ^ can be derived along the same lines as jSO] (Section 4, 
p. 364). It can be sketched as follows. First, it is clear that the transition probabil- 
ities of the proposed algorithm satisfy Equations 3.5-3.9 in 30j (p. 361-362). Next, 
in order to ensure the irreductibility of the Markov chain, the density of the process 
have to be hereditary (Definition 3.1 in |30], p. 360), which is detailed just below. 
Then by adapting the proof of Corollary 2 in Tierney (IHI], Section 3.1, p. 1713), it 
follows that the chain is Harris recurrent. 

In our context, the herediraty property is a direct consequence of the local stabil- 
ity of the process. Local stability means that the Energy variation can be controled 
under local modifications of a configuration. The following result can be used. 

Lemma 1 (local stability) Under the same conditions as in propositionUl there 
exists a constant C{Kj, A, d, D) so that 



where the inequality holds for all configurations and all Xn+i G with type r„_|_i. 
More specifically, the constant is equal to 



exp{-9H{ip)) ^ hnjip) 
Z ~ Z 



H{v) = J2 n x,)\J{t,, r,) + A 5^(|Dir(x,)| - A,fT{A,J, 




H{ip U (x„+i, r„+i)) - H{ip) I < C{K.j, A, rf, D) 



iV(A^ + l) D 
2 sin(/?) 
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D 
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C{K.j,X.d,D) 



2 sin(/?) 



where N = 2'k / (3 the maximum number of neighbors of any cell, and (3 is the mini- 
mum angle of all Delaunay's triangles. 

The proof of Lemma d is given in the Appendix. The restriction to J bounded 
and the other minor assumptions are necessary for technical reasons, and they may 
appear unreahstic as biological patterns are concerned. Nevertheless the situations 
encountered in the remainder of the article will always check these conditions. 

4 Results 

4.1 Simulation of biological patterns 

In this section, we report simulation results obtained with a finite set of marks 
M = {ti,T2,te}. As the GG model does, we show that the model has the abil- 
ity to reproduce at least three kinds of biological relevant patterns: Checkerboard, 
Cell Sorting and Engulfment. Checkerboard pattern formation was investigated in 
a simulation study of Honda et al. |32i about the sexual maturation of the avian 
oviduct epithelium. Cell sorting is a standard pattern of mixed heterotypic aggre- 
gates. Experimental observations of this phenomena were reported by Takeuchi et 
al. and Armstrong [Tj. Engulfment of a tissue by another one was studied by 
Armstrong and Foty et al. This phenomenon is a direct consequence of 

adhesion processes between the two cell types and the extracellular medium. 

The two marks ti and T2 represent "active cell types" with distinct phenotypes 
responsible for the adhesion process. In addition, active cells are surrounded by an 
extracellular medium modeled by cells of type te- These three types are similar to 
the /, d and M types of Glazier and Graner [22] • Simulations were generated from 
the Metropolis algorithm presented in the previous section. A unique configuration 
was used to initialize all the simulations. It consisted of about 1000 randomly located 
active cells, and the configuration displayed in Figure ^ In this configuration, the 
marks were also random. The target areas for active cells were equal to A^-^ = 
Ar^ = 5.10"^. Under equilibrium, configurations were expected to consist of about 
7r/5.10"'^ ~ 628 cells. No area constraint affected the te cells and we set Arj^ = 
— 1. The adhesion strength parameter 6 was fixed to ^ = 10, and we adjusted 
the interaction intensities so that the surface tensions were of same order as those 
considered by Granier and Glazier [22] . 

Checkerboard patterns can be interpreted as arising from negative surface ten- 
sions. In the GG model. Checkerboard patterns were generated using parameter 
settings that corresponded to surface tensions around 712 = —3. Figure [21 displays 



the configuration obtained after 50,000 cycles of tlie Metropolis-Hastings algoritlim, 
wliere tfie interaction intensities were fixed at J(ri, T2) = 0, J(ri, ri) = J{t2, T2) = 1 
and J{te, ti) = J{te, T2) = 0. Tliese interaction intensities corresponded to a sur- 
face tension equal to 712 = — 1 which was of the same order as the one used in the 
GG model. 

In contrast, Cell Sorting patterns arise from positive surface tensions between 
active cells. In the GG model. Cell Sorting patterns were generated using parameter 
settings that corresponded to surface tensions around 712 = -l-3. In our model, 
simulations were conducted using the following interaction intensities: J(ti, T2) = 1, 
J{ti,ti) = J{t~2,t~2) = and J{te,ti) = J{te,T2) = 0. This corresponded to the 
value 712 = +1. The configuration obtained after 50,000 steps is displayed in Figure 

m 

Simulations of Engulfment were conducted using the following parameters: 
J{n,r2) = 1, J{ri,Ti) = J{t2,T2) = 0, J{te,ti) = 0, J{te,T2) = 1. These 
interaction intensities provided positive surface tensions between active cells, which 
then tended to form clusters. The fact that J{te, T2) was greater than J{te, ti) 
ensured that ti cells were more likely to be close to the extracellular medium and 
to surround the T2 cells. The results are displayed in Figure El 

4.2 Statistical estimation of the adhesion strength parame- 
ter 

In this section, we study the influence of varying the adhesion strength parameter 6 
on simulation results, and we summarize the performances of the maximum pseudo- 
likelihood estimator 6. 

To assess the influence of 6 on pattern simulations, three values were tested: 
^ = 1, ^ = 5 and 6 = 10. The results are presented for simulations of Checkerboard 
and Cell Sorting patterns. In both cases, the interaction intensities were setting as 
in section 03 and provided 712 = —1 in Checkerboard simulations and 712 = +1 in 
Cell Sorting simulations. 

In all cases we ran the Metropolis algorithm for 50,000 steps, which were enough 
to provide a fiat profile of Energy variation. The three final configurations, in both 
Checkerboard and Cell Sorting, are displayed in Figure Either for Checkerboard 
or for Cell Sorting simulations, we can see a gradual evolution as 6 increases. For 
6 = 1, the marks are randomly distributed, for = 5 a small inhibition is visible in 
the Checkercoard simulation while small clusters appear in the Cell Sorting pattern. 
Finally, for 6 = 10 the stronger inhibition between cells with the same types provides 
a more pronounced checkerboard pattern, and larger clusters are obtained in Cell 
Sorting. 



We studied the statistical performances of 6 through simulations of Checkerboard 
and Cell Sorting patterns. The interaction intensities were the same as previously. 
For each value of 6, 50 replicates were generated from which the mean and the 
variance of 6 was estimated. Table El summarizes the results obtained for 6 in 
the range (1 . . .20). For Cell Sorting, the bias was weak for all values of 6, while 
for Checkerboard the bias seemed to be slightly higher. The results were similar 
regarding the variance. It was higher for Checkerboard than for Cell Sorting. In 
addition the variance increased as 9 increased. 

4.3 Real data 

Estimation of the adhesion strength was also performed on real data. We used 
survivin and beta-catenin markers in the context of medulloblastoma jHSj. These 
markers are known to be implicated in complexes that regulate adhesion between 
contiguous cells. An image analysis, analogous to the analysis performed in jHEj, has 
been achieved to extract the location of each cell nuclei and the level of expression 
of markers in each cell. The modeled tissue is displayed in Figure IHl 

We considered that the data were relevant to a Cell Sorting pattern, and the cor- 
responding interaction parameters were then used. The estimate of 6 was computed 
as ^ ~ 4.79 which fall into the range of well-fitted values. 

5 Discussion 

In this paper, a model based on marked point processes theory and inspired from 
[in] has been studied. The new model proposes a continuous extension of the GG 
model by considering continuous instead of discrete cells, and it overcomes some 
limitations of the discrete model. First, in the discrete simulation algorithm cells are 
not constrained to be simply connected (i.e., they may be divided in non-connected 
components). Graner and Glazier solved this problem by choosing appropriate initial 
configurations and small temperatures. Another issue related to the discrete model is 
that the discretization scale influences the choice of the parameters A and T which 
in turn contribute to the acceptance/rejection probabilities. Finally the discrete 
algorithm was lacking good convergence properties. More specifically, the Markov 
chain generated by the discrete algorithm might not be ergodic and might also be 
influenced by the discretization scale. Again, these drawbacks did not impact the 
original results, because the authors chose their initial configurations so that the 
desired output was produced. 

The algorithm proposed in section El provides convergent simulations whatever 
the initial configuration and is also considerably faster. In this study, we have shown 



that the new model was able to reproduce biologically relevant cell patterns such as 
Checkerboard, Cell Sorting and Engulfment. Furthermore, the model has been built 
so that it include the strength of cell-cell adhesion. We have proposed and validated 
an inference procedure based on the pseudo-likelihood. The statistical errors were 
low in Cell Sorting simulations. In Checkerboard simulations, bias and variance 
were slighlty higher than for Cell Sorting but still reasonable in the range of tested 
parameters. 

Further improvements of this approach would require a deeper mathematical 
study which is beyond the scope of this study. In particular, the theory of marked 
point processes makes it possible to establish theoretical consistency results for 6. 
Furthermore, the other interaction parameters can also be estimated in the same 
way as was. Although we did not assessed the performances of these estimators, 
we believe that they would be useful for analyzing tissues arrays, as generated by 
high-throughtput cancer studies [HTj . 
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Appendix 



In this section, we give a proof of Lemma ^ This proof needs some geometrical 
results that are summarized in the following Lemma. 

Lemma 2 Let ip = {(xi, ri), . . . , r„)} be a configuration. Assume that there 
exists {d, D) such as < d < D < +00 and 



Xi ^ Xj ^ d < \xi ~ Xj\ < D. 



Then we have: 



1. the minimum angle in all Delaunay's triangles of (p is greater than 

d' 



j3 = arccos 1 



2. n{xi) < where n{xi) is the number of neighbors of Xi in the configuration 

3. iDirfx,- n < 



D 



on — sm(/3)' 

2 



4. |Dir(x.)|<7r(^ 



Proof of Lemma [2] 

1. This result is derived from the law of cosines 

2. comes from item 1. and the fact that each Delaunay's neighbor belongs to 
two Delaunay's triangles. 

3. The radius of a circumcircle to a triangle is equal to a/ sin(Q;) where a is the 
length of an edge an a is the opposite angle to the edge. We deduce that the radius 
of all circumcircle is less than D/(2 sin(/3)). Then we have for all z ~ j 

D 

Dir nxj- < —7^ 
sm(p) 

4. From 3., we deduce that Dir(xj) is included in the circle of center Xi and of 
radius D/ sin(/5). 



Proof of Lemma [T] 

This proof can be derived along the same lines as [2^1 • Let = {(xi, Xj), . . . , (x„, r„)} 
be a (marked) configuration and (x„+i,r„+i) be a (marked) cell. We write 



V?' = V? U (x„+i,r„+i). 



As Dirichlet cells are associated with a given configuration, we write Dir(xj) 
(resp. Dir'(xj)) the Dirichlet cell associated with the cell Xi in the configuration ip 
(resp. ip'). Analogously, Dir^Xinxj) (resp. Dir'(xi fl x-,)) characterizes the Dirichlet 
edge in the configuration ip (resp. ip'). 

Then, according to the definition of H in Equation 01 the difference between 
H{ip') and H{ip) is given by 

H{ip')-H{^) = 5];|Dir'(x.nx,)|J(r„r,) + Aj](|Dir'(x.)|-A.j2r(A.J 

J2 |Dir(x. n x,)\JiT„ r,) - A 5^(|Dir(a:.)| - A^J'ViA^J 



Using geometrical properties of Dirichlet tesselation, the previous expression can 
be simplified in a sum of four terms. 



+A(|Dir'|(x„+i)-A.„^j2r(A.„_,J 

+ {\DiT'{xinxj) \ - \Bii{xir]Xj)\)J{ri,Tj) 



Triangles(i,j,n+1) 



+A Y mAxi)\ - AJ^r(AJ - (|Dir(x,)| - A^fT{A^^] 



i;i~n+l 

The first term corresponds to the interaction of the additional cell Xn+i with its 
neighbors. The second term is explained by the shape constraint for the additional 
cell Xn+i- The third term stand for the difference between the interactions in ip' 
and in ip. The sum runs over all triangles n+1) that belong to the Delaunay's 
configuration. Using geometrical properties about the insertion of a Dirichlet cell, 
only interactions between two cells that are both in the neighborhood of X'h^i are 
non-zero. The fourth term represents the difference between shape constraints in ip' 
and in ip. According to geometrical properties, the sum runs over all cells in the 
neighborhood of Xn+i- 

Each term can be controled independently. First, from LemmaElitem 3, we have 
|Dir'(a;j fl < D/ sin(/5). Furthermore, Lemma|21item 2 gives that the number 

of neighbors of Xn+i is less than A^, which leads to the following result 



|Dir'(a;i n x„+i)|J(ri,r„+i) 
Next, directly from Lemma |21 item 4, we have 



< N-^Kj. (5) 



|A(|Dir'|(x„+0 - A.„,J^r(A,,,j| < Xn' (^^^^ 



(6) 



From Lemma 121 item 3, we have |Dir'(xj fl x„+i)| < D / sin(/5). In addition, from 
Lemma El item 2, we extract that the number of pairs in the neighborhood of Xn+i 
is less than A^(A^ — l)/2, which leads to the following result 



^ (|Dir'(xi nxj)| - \Dii{xinxj)\)J{n,Tj] 

Triangles(i,j,n+1) 

Finally, from Lemma |21 item 2 and |21 item 4 we have 



^ N{N - 1) D 



A (|Dir'(a;0|-A.J'r(A.J-(|Dir(xO|-A.jm 



Combining Equations EHHl we obtain that 



D 



2 sin(/3) 



(8) 



sin(/?)' 



2 sin(/5) 







Checkerboarder 


Cell Sorting 






Mean 


Variance 


Mean 


Variance 




1 


0.98 


0.70 


1.03 


0.4 




3 


3.14 


0.66 


3.01 


0.51 




: 5 


5.01 


0.57 


4.94 


0.94 




: 8 


8.20 


1.07 


8.01 


0.81 




10 


10.47 


1.20 


9.80 


1.00 




12 


12.28 


1.81 


12.05 


1.09 




15 


14.58 


2.22 


15.03 


1.20 




20 


20.44 


3.55 


20.08 


2.98 



Table 1: Mean and Variance of 9, maximum of Pseudo-likelihood, for Checkerboard 
and Cell Sorting simulations. The evaluation was achievied using 50 simulations of 
each case. 
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Figure 1: The initial configuration for simulating Checkerboard, Cell Sorting and En- 
gulfment patterns. It consists of about 1000 active cells surrounded by an extracellular 
medium. The active cells are randomly located in the unit sphere, and their types are 
randomly sampled from M. Cells of type ri are colored grey while cells of type T2 are 
colored in black. 



-1.0 -0.5 0.0 05 1.0 10000 20000 30000 40000 50000 

X 1:50001 

(a) (b) 

Figure 2: Checkerboard simulation. The interaction intensities were chosen as follows: 
J{ri,Ti) = 1, J{t2,T2) = 1, J{ti,T2) = 0, J(ti,te) = and J{t2,te) = 0. (a) The 
configuration obtained after 50,000 iterations, (b) The decrease of the Energy as a function 
of the iteration steps. 




(a) (b) 

Figure 3: Cell Sorting simulation. The interaction intensities were chosen as follows: 
J(ri,ri) = 0, J{t2,T2) = 0, J(ri,r2) = 1, J{ti,te) = and J{t2,te) = 0. (a) The 
configuration obtained after 50,000 iterations, (b) The decrease of the Energy as a function 
of the iteration steps. 




Figure 4: Engulfment simulation. The interaction intensities were chosen as follows: 
J{ri,Ti) = 0, J{t2,T2) = 0, J{ti,T2) = 1, J{ti,te) = and J{t2,te) = 1. (a) The 
configuration obtained after 50,000 iterations, (b) The decrease of the Energy as a function 
of the iteration steps. 
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Figure 5: Influence of 9 in Checkerboard simulations. Final configurations using three 
different values for 9. Simulations gradually corresponds to either a Checkerboard or large 
clusters. 




Figure 6: Real data. The statistical procedure was conducted using the interactions 
parameters for cell sorting patterns: J(ri,ri) = 1, J(t2, T2) = 1, J{ti,T2) — 0, 
J{ti,te) = and J(t2, te) = 0. We obtained 9 ^ 4.79. 



